
import aggregateGeneCoverage
import geneSummary
import sys

def printUsage(args, usage, minArgs=2):
  if len(args) < minArgs or args[1].startswith("-h") or args[1].startswith("--h"):
    print usage
    sys.exit(1)
    
usageString = """
Usage: python %s SUMMARY_FILE(S) OUTPUT_NAME

where SUMMARY_FILE(S) is a comma-separated list of
mapping summary files, and output will be written
to OUTPUT_NAME.exon-intron-coverage.txt and
OUTPUT_NAME.normalized-gene-summary.txt.
"""

if __name__ == "__main__":
    printUsage(sys.argv, usageString % sys.argv[0], 3)
    summaryFiles = sys.argv[1].split(",")
    outputName = sys.argv[2]
    exonIntronFile = outputName + ".exon-intron-coverage.txt"
    geneSummaryFile = outputName + ".normalized-gene-summary.txt"

    print >> sys.stderr, "# Summary files: ", summaryFiles
    print >> sys.stderr, "# Exon/intron coverage output: ", exonIntronFile
    print >> sys.stderr, "# Gene expression output: ", geneSummaryFile
    
    exonIntronStream = open(exonIntronFile, "w")
    humanMappedReads = aggregateGeneCoverage.aggregateAllChromosomes(summaryFiles, aggregateGeneCoverage.KNOWN_GENE_TABLE, exonIntronStream)
    exonIntronStream.close()

    stdout = sys.stdout
    sys.stdout = open(geneSummaryFile, "w")
    if sum(humanMappedReads) > 0:
      print >> sys.stderr, "# Counted %d mapped reads; calling gene summary with norm = %.7f" % (sum(humanMappedReads), 1000000.0 / sum(humanMappedReads))
      geneSummary.geneSummary(exonIntronFile, 1000000.0 / sum(humanMappedReads))
    else:
      print >> sys.stderr, "# Counted 0 mapped reads; not calling gene summary"
    sys.stdout.close()
    sys.stdout = stdout

